load nursery.mat;
reset(RandStream.getDefaultStream)
nurseryTest = nursery(randperm(size(nursery,1)),:);

% Partition training data
train = nurseryTest(1:size(nursery,1)/2,:);

% Partition test data
test = nurseryTest(size(nursery,1)/2+1:end,:);

numXValues = size(nurseryTest,2) -1;

%find P(Y)
yFreq = zeros(5,1);
for trainIndex=1:5
    yFreq(trainIndex,1) = sum(train(:,numXValues+1)==trainIndex);
end
yProb = yFreq./size(train,1);

%trainResults has number of x BY possible values of x BY possible values of Y
%trainResults contain the conditional probabilities for each x (i.e., P(x = 1 | y = 2))
trainResults = zeros(numXValues,5,5);
for xIndex=1:numXValues
    for xvalue=1:5
        for yvalue=1:5
            trainResults(xIndex,xvalue,yvalue) = sum(train(:,numXValues+1)==yvalue&train(:,xIndex)==xvalue) /sum(train(:,numXValues+1)==yvalue);
        end
    end
end

%get probability distribution over Y given x1:x8
testProbabilities =  ones(size(test,1),5);
for testIndex = 1:size(test,1)
    for yIndex = 1:5
        for xIndex = 1:numXValues        
            testProbabilities(testIndex,yIndex) = testProbabilities(testIndex,yIndex) * trainResults(xIndex,test(testIndex,xIndex),yIndex);
        end
        testProbabilities(testIndex,yIndex) = testProbabilities(testIndex,yIndex) * yProb(yIndex,1);
    end
end

%classified Y = greatest probability of Y given x1:x8
[~,testIndexResults] = max(testProbabilities');
testIndexResults = testIndexResults';

%calculate accuracy
A = [testIndexResults test(:,numXValues+1)];
accuracy = sum(A(:,1)==A(:,2)) / size(A,1);

%calculate log likelihood of training data
logLikelihood = 0;
for index=1:size(train,1)
    for dimensionIndex=1:numXValues
        logLikelihood = logLikelihood + log(trainResults(dimensionIndex,train(index, dimensionIndex),train(index, numXValues + 1)));
    end
    logLikelihood = logLikelihood + log(yProb(train(index, numXValues + 1)));
end


